Setup R notebook Environment

knitr::opts_chunk$set(echo = TRUE)
Sys.setlocale("LC_CTYPE", "en_US.UTF-8")
[1] ""

Package load for Project

packages = c("sp","sf","tidyverse","tmap","rjson", "rgdal", "geojsonio", "GISTools",'spatialEco', "raster","dplyr") 
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p) 
  } 
  library(p,character.only = T) 
}
Loading required package: sp
Loading required package: sf
Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
Loading required package: tidyverse
-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
v ggplot2 3.1.0     v purrr   0.2.5
v tibble  1.4.2     v dplyr   0.7.8
v tidyr   0.8.2     v stringr 1.3.1
v readr   1.3.1     v forcats 0.3.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
Loading required package: tmap
Loading required package: rjson
Loading required package: rgdal
rgdal: version: 1.3-6, (SVN revision 773)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/terence.tan.2015/Documents/R/win-library/3.5/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/terence.tan.2015/Documents/R/win-library/3.5/rgdal/proj
 Linking to sp version: 1.3-1 
Loading required package: geojsonio

Attaching package: <U+393C><U+3E31>geojsonio<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:base<U+393C><U+3E32>:

    pretty

Loading required package: GISTools
Loading required package: maptools
Checking rgeos availability: TRUE
Loading required package: RColorBrewer
Loading required package: MASS

Attaching package: <U+393C><U+3E31>MASS<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:dplyr<U+393C><U+3E32>:

    select

Loading required package: rgeos
rgeos version: 0.4-2, (SVN revision 581)
 GEOS runtime version: 3.6.1-CAPI-1.10.1 
 Linking to sp version: 1.3-1 
 Polygon checking: TRUE 

Loading required package: spatialEco
Loading required package: raster

Attaching package: <U+393C><U+3E31>raster<U+393C><U+3E32>

The following objects are masked from <U+393C><U+3E31>package:MASS<U+393C><U+3E32>:

    area, select

The following object is masked from <U+393C><U+3E31>package:dplyr<U+393C><U+3E32>:

    select

The following object is masked from <U+393C><U+3E31>package:tidyr<U+393C><U+3E32>:

    extract

Importing Taiwan GDAM Subzone

Shapefile Object with Spatial data Import data for past 12 years

Requires translation of chinese character independetly (https://stat.ethz.ch/pipermail/r-sig-geo/2009-March/005323.html)

Taiwan follows EPSG 3826 (http://spatialreference.org/ref/epsg/twd97-tm2-zone-121/)

taiwan_ts_map_sp <- readOGR(dsn = "data/TAIWAN_TOWNSHIP", layer = "TOWN_MOI_1071226", stringsAsFactors=TRUE)
OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\terence.tan.2015\Desktop\dangy\data\TAIWAN_TOWNSHIP", layer: "TOWN_MOI_1071226"
with 368 features
It has 7 fields
#taiwan <- readShapePoly("data/taiwan_data/COUNTY_MOI_1070516.shp")
#taiwan@proj4string<- CRS( "+init=epsg:3826 +proj=longlat +ellps=WGS84 +no_defs")
#taiwan.union <- aggregate(taiwan)
df_dengue <- jsonlite::fromJSON("data/dengue_case.json") %>% 
  filter(as.Date(Onset_day) >= "2016-01-01" & as.Date(Onset_day)< "2017-01-01")
# Check for NA values for coordinates
 sum(is.na(df_dengue$Minimum_statistical_area_center_point_X))
[1] 0
 sum(is.na(df_dengue$Minimum_statistical_area_center_point_Y))
[1] 0
 df_dengue<- df_dengue[!is.na(df_dengue$Minimum_statistical_area_center_point_X),]
 df_dengue<- df_dengue[!(df_dengue$Minimum_statistical_area_center_point_X == 'None'),]
 
 # Type conversion
 df_dengue[, c(10,11,19,23,24)] <- sapply(df_dengue[, c(10,11,19,23,24)], as.numeric)
 
  # Transform into SF object
 sf_dengue <- st_as_sf(df_dengue, 
                       coords = c("Minimum_statistical_area_center_point_X",
                                  "Minimum_statistical_area_center_point_Y"),
                        crs = "+init=epsg:3826 +proj=longlat +ellps=WGS84 +no_defs",na.fail=FALSE)
sf_dengue <- na.omit(sf_dengue)
sf_dengue <- as(sf_dengue, 'Spatial')
sp_dengue <- as(sp_dengue, 'SpatialPoints')
tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(taiwan_ts_map_sp)+
       tm_borders(col = "grey40",alpha=0.5)+
  # prevent zoom too much
  tm_view(alpha = 1, set.zoom.limits = c(7,21))+
  tm_layout(basemaps = c('OpenStreetMap'))+
  tm_shape(sp_dengue)+
  tm_dots(col ="red",size = 0.01) 
As of version 2.0, basemaps and basemaps.alpha have to be called from tm_view

tmap_mode("plot")
tmap mode set to plotting

Converting the generic sp format into spatstat’s ppp format

dengue_ppp <- as(sp_dengue, "ppp")
summary(dengue_ppp) 
Planar point pattern:  320 points
Average intensity 70.15672 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 6 decimal places

Window: rectangle = [120.14635, 121.79824] x [22.411063, 25.172282] units
Window area = 4.56122 square units

Combining GP clinics points and the study area

tw_owin <- as(taiwan_ts_map_sp, "owin") 
plot(tw_owin) 

summary(tw_owin)
Window: polygonal boundary
1007 separate polygons (no holes)
enclosing rectangle: [118.13797, 124.56115] x [21.8956, 26.385278] units
Window area = 3.26423 square units
Fraction of frame area: 0.113
dengueTW_ppp = dengue_ppp[tw_owin]
summary(dengueTW_ppp)
Planar point pattern:  320 points
Average intensity 98.03217 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 6 decimal places

Window: polygonal boundary
1007 separate polygons (no holes)
enclosing rectangle: [118.13797, 124.56115] x [21.8956, 26.385278] units
Window area = 3.26423 square units
Fraction of frame area: 0.113
plot(dengueTW_ppp)

qt <- quadrat.test(dengueTW_ppp,                     
                   nx = 20, ny = 15)
Some expected counts are small; chi^2 approximation may be inaccurate
plot(dengueTW_ppp) 
plot(qt, add = TRUE, cex =.1) 

qt

Kernel Density Estimation

#dengueTW_ppp <- rescale(dengueTW_ppp, 1000, "km") 
kde_taiwan_bw <- density(dengueTW_ppp, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_taiwan_bw)

LS0tDQp0aXRsZTogImZlYXR1cmUta2RhIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KU2V0dXAgUiBub3RlYm9vayBFbnZpcm9ubWVudA0KYGBge3Igc2V0dXAsIGluY2x1ZGU9VFJVRSwgZXZhbD1UUlVFLG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUpDQpTeXMuc2V0bG9jYWxlKCJMQ19DVFlQRSIsICJlbl9VUy5VVEYtOCIpDQpgYGANCg0KI1BhY2thZ2UgbG9hZCBmb3IgUHJvamVjdA0KYGBge3J9DQpwYWNrYWdlcyA9IGMoInNwIiwic2YiLCJ0aWR5dmVyc2UiLCJ0bWFwIiwicmpzb24iLCAicmdkYWwiLCAiZ2VvanNvbmlvIiwgIkdJU1Rvb2xzIiwnc3BhdGlhbEVjbycsICJyYXN0ZXIiLCJkcGx5ciIsICJzcGF0c3RhdCIsDQogICAgICAgICAgICAgJ21hcHRvb2xzJywgJ2xlYWZsZXQnKSANCmZvciAocCBpbiBwYWNrYWdlcyl7DQogIGlmKCFyZXF1aXJlKHAsIGNoYXJhY3Rlci5vbmx5ID0gVCkpew0KICAgIGluc3RhbGwucGFja2FnZXMocCkgDQogIH0gDQogIGxpYnJhcnkocCxjaGFyYWN0ZXIub25seSA9IFQpIA0KfQ0KYGBgDQoNCiNJbXBvcnRpbmcgVGFpd2FuIEdEQU0gU3Viem9uZQ0KU2hhcGVmaWxlIE9iamVjdCB3aXRoIFNwYXRpYWwgZGF0YQ0KSW1wb3J0IGRhdGEgZm9yIHBhc3QgMTIgeWVhcnMNCg0KUmVxdWlyZXMgdHJhbnNsYXRpb24gb2YgY2hpbmVzZSBjaGFyYWN0ZXIgaW5kZXBlbmRldGx5IChodHRwczovL3N0YXQuZXRoei5jaC9waXBlcm1haWwvci1zaWctZ2VvLzIwMDktTWFyY2gvMDA1MzIzLmh0bWwpDQoNClRhaXdhbiBmb2xsb3dzIEVQU0cgMzgyNiAoaHR0cDovL3NwYXRpYWxyZWZlcmVuY2Uub3JnL3JlZi9lcHNnL3R3ZDk3LXRtMi16b25lLTEyMS8pDQpgYGB7cn0NCnRhaXdhbl90c19tYXBfc3AgPC0gcmVhZE9HUihkc24gPSAiZGF0YS9UQUlXQU5fVE9XTlNISVAiLCBsYXllciA9ICJUT1dOX01PSV8xMDcxMjI2Iiwgc3RyaW5nc0FzRmFjdG9ycz1UUlVFKQ0KDQojdGFpd2FuIDwtIHJlYWRTaGFwZVBvbHkoImRhdGEvdGFpd2FuX2RhdGEvQ09VTlRZX01PSV8xMDcwNTE2LnNocCIpDQojdGFpd2FuQHByb2o0c3RyaW5nPC0gQ1JTKCAiK2luaXQ9ZXBzZzozODI2ICtwcm9qPWxvbmdsYXQgK2VsbHBzPVdHUzg0ICtub19kZWZzIikNCiN0YWl3YW4udW5pb24gPC0gYWdncmVnYXRlKHRhaXdhbikNCg0KDQpkZl9kZW5ndWUgPC0ganNvbmxpdGU6OmZyb21KU09OKCJkYXRhL2Rlbmd1ZV9jYXNlLmpzb24iKSAlPiUgDQogIGZpbHRlcihhcy5EYXRlKE9uc2V0X2RheSkgPj0gIjE5OTgtMDEtMDEiICYgYXMuRGF0ZShPbnNldF9kYXkpPCAiMTk5OS0wMS0wMSIpDQoNCiMgQ2hlY2sgZm9yIE5BIHZhbHVlcyBmb3IgY29vcmRpbmF0ZXMNCiBzdW0oaXMubmEoZGZfZGVuZ3VlJE1pbmltdW1fc3RhdGlzdGljYWxfYXJlYV9jZW50ZXJfcG9pbnRfWCkpDQogc3VtKGlzLm5hKGRmX2Rlbmd1ZSRNaW5pbXVtX3N0YXRpc3RpY2FsX2FyZWFfY2VudGVyX3BvaW50X1kpKQ0KIGRmX2Rlbmd1ZTwtIGRmX2Rlbmd1ZVshaXMubmEoZGZfZGVuZ3VlJE1pbmltdW1fc3RhdGlzdGljYWxfYXJlYV9jZW50ZXJfcG9pbnRfWCksXQ0KIGRmX2Rlbmd1ZTwtIGRmX2Rlbmd1ZVshKGRmX2Rlbmd1ZSRNaW5pbXVtX3N0YXRpc3RpY2FsX2FyZWFfY2VudGVyX3BvaW50X1ggPT0gJ05vbmUnKSxdDQogDQogIyBUeXBlIGNvbnZlcnNpb24NCiBkZl9kZW5ndWVbLCBjKDEwLDExLDE5LDIzLDI0KV0gPC0gc2FwcGx5KGRmX2Rlbmd1ZVssIGMoMTAsMTEsMTksMjMsMjQpXSwgYXMubnVtZXJpYykNCiANCiAgIyBUcmFuc2Zvcm0gaW50byBTRiBvYmplY3QNCiBzZl9kZW5ndWUgPC0gc3RfYXNfc2YoZGZfZGVuZ3VlLCANCiAgICAgICAgICAgICAgICAgICAgICAgY29vcmRzID0gYygiTWluaW11bV9zdGF0aXN0aWNhbF9hcmVhX2NlbnRlcl9wb2ludF9YIiwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAiTWluaW11bV9zdGF0aXN0aWNhbF9hcmVhX2NlbnRlcl9wb2ludF9ZIiksDQogICAgICAgICAgICAgICAgICAgICAgICBjcnMgPSAiK2luaXQ9ZXBzZzozODI2ICtwcm9qPWxvbmdsYXQgK2VsbHBzPVdHUzg0ICtub19kZWZzIixuYS5mYWlsPUZBTFNFKQ0Kc2ZfZGVuZ3VlIDwtIG5hLm9taXQoc2ZfZGVuZ3VlKQ0Kc3BfZGVuZ3VlIDwtIGFzKHNmX2Rlbmd1ZSwgJ1NwYXRpYWwnKQ0KYGBgDQoNCmBgYHtyfQ0Kc3BfZGVuZ3VlIDwtIGFzKHNwX2Rlbmd1ZSwgJ1NwYXRpYWxQb2ludHMnKQ0KYGBgDQoNCmBgYHtyfQ0KdG1hcF9tb2RlKCJ2aWV3IikNCnRtX3NoYXBlKHRhaXdhbl90c19tYXBfc3ApKw0KICAgICAgIHRtX2JvcmRlcnMoY29sID0gImdyZXk0MCIsYWxwaGE9MC41KSsNCiAgIyBwcmV2ZW50IHpvb20gdG9vIG11Y2gNCiAgdG1fdmlldyhhbHBoYSA9IDEsIHNldC56b29tLmxpbWl0cyA9IGMoNywyMSkpKw0KICB0bV9sYXlvdXQoYmFzZW1hcHMgPSBjKCdPcGVuU3RyZWV0TWFwJykpKw0KICB0bV9zaGFwZShzcF9kZW5ndWUpKw0KICB0bV9kb3RzKGNvbCA9InJlZCIsc2l6ZSA9IDAuMDEpIA0KDQp0bWFwX21vZGUoInBsb3QiKQ0KYGBgDQoNCg0KIyBDb252ZXJ0aW5nIHRoZSBnZW5lcmljIHNwIGZvcm1hdCBpbnRvIHNwYXRzdGF0J3MgcHBwIGZvcm1hdCANCg0KYGBge3IgZWNobz1UUlVFLCBldmFsPVRSVUV9DQpkZW5ndWVfcHBwIDwtIGFzKHNwX2Rlbmd1ZSwgInBwcCIpDQpzdW1tYXJ5KGRlbmd1ZV9wcHApIA0KYGBgDQoNCiMgQ29tYmluaW5nICBHUCBjbGluaWNzIHBvaW50cyBhbmQgdGhlIHN0dWR5IGFyZWEgDQpgYGB7ciBlY2hvPVRSVUUsIGV2YWw9VFJVRX0NCnR3X293aW4gPC0gYXModGFpd2FuX3RzX21hcF9zcCwgIm93aW4iKSANCnBsb3QodHdfb3dpbikgDQpzdW1tYXJ5KHR3X293aW4pDQoNCmRlbmd1ZVRXX3BwcCA9IGRlbmd1ZV9wcHBbdHdfb3dpbl0NCnN1bW1hcnkoZGVuZ3VlVFdfcHBwKQ0KcGxvdChkZW5ndWVUV19wcHApDQpgYGANCmBgYHtyfQ0KcXQgPC0gcXVhZHJhdC50ZXN0KGRlbmd1ZVRXX3BwcCwgICAgICAgICAgICAgICAgICAgICANCiAgICAgICAgICAgICAgICAgICBueCA9IDIwLCBueSA9IDE1KQ0KYGBgDQoNCmBgYHtyfQ0KcGxvdChkZW5ndWVUV19wcHApIA0KcGxvdChxdCwgYWRkID0gVFJVRSwgY2V4ID0uMSkgDQpgYGANCg0KYGBge3J9DQpxdA0KYGBgDQoNCg0KIyBLZXJuZWwgRGVuc2l0eSBFc3RpbWF0aW9uDQoNCmBgYHtyIGVjaG89VFJVRSwgZXZhbD1UUlVFfQ0KZGVuZ3VlVFdfcHBwIDwtIHJlc2NhbGUoZGVuZ3VlVFdfcHBwLCAxMDAwLCAia20iKSANCiNrZGVfdGFpd2FuX2J3IDwtIGRlbnNpdHkoZGVuZ3VlVFdfcHBwLCBzaWdtYT1idy5kaWdnbGUsIGVkZ2U9VFJVRSwga2VybmVsPSJnYXVzc2lhbiIpDQoNCiNwbG90KGtkZV90YWl3YW5fYncpDQpgYGANCg0KYGBge3J9DQpwbG90KGtkZV90YWl3YW5fYncpDQpgYGANCg0K